sun <- read_rds("~/GoogleDrive/PhD/Paper Three/mob_multicity.rds") %>% filter(codigo_destino != codigo_residencia & dia == "Sunday") %>% na.omit()
wed <- read_rds("~/GoogleDrive/PhD/Paper Three/mob_multicity.rds") %>% filter(codigo_destino != codigo_residencia & dia == "Wednesday") %>% na.omit()

Part 1: Basic Gravity Model - Density, Distance, City-as-Level

s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = sun)
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = wed)

Part 2: Can more complex models wtih land cover and/or weather data improve on the best basic model?

Land Cover Model 1 - Urban Fabric Continuity Only

s1_lc <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + scale(uf_cont_destino) + scale(uf_cont_residencia), data = sun)
w1_lc <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + scale(uf_cont_destino) + scale(uf_cont_residencia), data = wed)

Land Cover Model 2 - All Land Cover Vars

lc2_formula <- as.formula(paste0("log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + ", paste0("scale(",names(sun)[12:23],")", collapse = " + ")))

s2_lc <- lmer(formula = lc2_formula, data = sun)
w2_lc <- lmer(formula = lc2_formula, data = wed)

Weather Model

s3_wea <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip, data = sun)
w3_wea <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip, data = wed)

Land Cover & Weather Model

all_formula <- as.formula(paste0("log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip + ", paste0("scale(",names(sun)[12:23],")", collapse = " + ")))

s4_all <- lmer(formula = all_formula, data = sun)
w4_all <- lmer(formula = all_formula, data = wed)

Comparison of Complex Models to Basic Model

AIC(s_basic,s1_lc,s2_lc,s3_wea,s4_all)
##         df      AIC
## s_basic  7 552433.7
## s1_lc    9 552144.1
## s2_lc   19 537242.1
## s3_wea  10 551799.6
## s4_all  22 536568.5
AIC(w_basic,w1_lc,w2_lc,w3_wea,w4_all)
##         df     AIC
## w_basic  7 1144243
## w1_lc    9 1143452
## w2_lc   19 1129340
## w3_wea  10 1142632
## w4_all  22 1127414

Model with All Land Cover and Weather Data Performs Best

summary(s4_all)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) +  
##     I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip +  
##     scale(uf_cont_destino) + scale(parks_destino) + scale(agriculture_destino) +  
##     scale(natural_destino) + scale(industry_destino) + scale(water_destino) +  
##     scale(uf_cont_residencia) + scale(parks_residencia) + scale(agriculture_residencia) +  
##     scale(natural_residencia) + scale(industry_residencia) +  
##     scale(water_residencia)
##    Data: sun
## 
## REML criterion at convergence: 536524.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2995 -0.6948 -0.0535  0.6546  4.3021 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  city     (Intercept) 0.1374   0.3707  
##  Residual             0.4590   0.6775  
## Number of obs: 260426, groups:  city, 6
## 
## Fixed effects:
##                                 Estimate Std. Error  t value
## (Intercept)                    2.616e+01  1.985e-01  131.821
## scale(dens_destino)           -1.873e-01  2.400e-03  -78.037
## scale(dens_residencia)         1.227e-02  2.379e-03    5.158
## log(dist)                     -4.536e+00  3.259e-02 -139.216
## I(log(dist)^2)                 2.215e-01  2.068e-03  107.117
## temp                           1.008e-03  8.335e-04    1.210
## I(temp^2)                     -1.162e-04  1.806e-05   -6.433
## precip                        -2.297e-04  2.346e-05   -9.790
## scale(uf_cont_destino)         2.054e-02  1.766e-03   11.629
## scale(parks_destino)           5.549e-03  1.655e-03    3.352
## scale(agriculture_destino)     9.500e-02  2.316e-03   41.024
## scale(natural_destino)         8.498e-02  2.467e-03   34.446
## scale(industry_destino)        1.683e-02  1.919e-03    8.766
## scale(water_destino)           9.044e-02  1.556e-03   58.132
## scale(uf_cont_residencia)     -1.456e-02  1.773e-03   -8.213
## scale(parks_residencia)        3.121e-02  1.646e-03   18.954
## scale(agriculture_residencia)  6.953e-02  2.150e-03   32.346
## scale(natural_residencia)      8.835e-02  2.233e-03   39.564
## scale(industry_residencia)     1.325e-02  1.854e-03    7.146
## scale(water_residencia)        2.749e-02  1.508e-03   18.229
## 
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
summary(w4_all)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) +  
##     I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip +  
##     scale(uf_cont_destino) + scale(parks_destino) + scale(agriculture_destino) +  
##     scale(natural_destino) + scale(industry_destino) + scale(water_destino) +  
##     scale(uf_cont_residencia) + scale(parks_residencia) + scale(agriculture_residencia) +  
##     scale(natural_residencia) + scale(industry_residencia) +  
##     scale(water_residencia)
##    Data: wed
## 
## REML criterion at convergence: 1127370
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6722 -0.7201 -0.0640  0.6798  4.3742 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  city     (Intercept) 0.2811   0.5302  
##  Residual             0.4435   0.6659  
## Number of obs: 556646, groups:  city, 6
## 
## Fixed effects:
##                                 Estimate Std. Error  t value
## (Intercept)                    2.622e+01  2.353e-01  111.429
## scale(dens_destino)           -2.262e-01  1.572e-03 -143.901
## scale(dens_residencia)         2.899e-02  1.596e-03   18.160
## log(dist)                     -4.681e+00  2.283e-02 -205.074
## I(log(dist)^2)                 2.415e-01  1.410e-03  171.240
## temp                           2.065e-02  6.044e-04   34.160
## I(temp^2)                     -5.218e-04  1.335e-05  -39.082
## precip                         9.593e-05  1.791e-05    5.355
## scale(uf_cont_destino)        -8.779e-04  1.218e-03   -0.720
## scale(parks_destino)          -7.387e-02  1.127e-03  -65.567
## scale(agriculture_destino)     3.981e-03  1.502e-03    2.650
## scale(natural_destino)         1.343e-02  1.577e-03    8.518
## scale(industry_destino)        3.213e-02  1.251e-03   25.672
## scale(water_destino)           3.356e-02  1.007e-03   33.321
## scale(uf_cont_residencia)     -1.010e-02  1.169e-03   -8.641
## scale(parks_residencia)        1.091e-02  1.129e-03    9.667
## scale(agriculture_residencia)  3.713e-02  1.447e-03   25.663
## scale(natural_residencia)      7.813e-02  1.514e-03   51.590
## scale(industry_residencia)     1.061e-02  1.231e-03    8.622
## scale(water_residencia)        1.562e-02  9.661e-04   16.163
## 
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
s_complex <- s4_all
w_complex <- w4_all

Part 3: Leave Out City, Make Predictions, Compare to Observation

s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(sun,city != "B"))
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(wed,city != "B"))
s_complex <- lmer(formula = all_formula, data = filter(sun,city != "B"))
w_complex <- lmer(formula = all_formula, data = filter(wed,city != "B"))
  
sun_pred <- sun %>% filter(city == "B") %>% mutate(flujo_basic = exp(predict(s_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(s_complex, ., allow.new.levels = T)))
wed_pred <- wed %>% filter(city == "B") %>% mutate(flujo_basic = exp(predict(w_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(w_complex, ., allow.new.levels = T)))

map <- suppressMessages(st_read("~/GoogleDrive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp")) %>% filter(ID_GRUPO %in% unique(sun_pred$codigo_destino))
## Reading layer `celdas_marzo_2020' from data source 
##   `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
sun_res <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)

sun_des <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)

sun_geo <- bind_rows(sun_res,sun_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)

wed_res <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)

wed_des <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)

wed_geo <- bind_rows(wed_res,wed_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)

p1 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = sun_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")

p2 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")

p3 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Complex Model")

p4 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = wed_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")

p5 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")

p6 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Complex Model")

Barcelona

Sundays

Wednesdays

Madrid

## Reading layer `celdas_marzo_2020' from data source 
##   `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator

Sundays

Wednesdays

Valencia

## Reading layer `celdas_marzo_2020' from data source 
##   `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator

Sundays

Wednesdays

Sevilla

## Reading layer `celdas_marzo_2020' from data source 
##   `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator

Sundays

Wednesdays